## This file is for analysing spatial panel data set "txData" using R package "plm" and "splm"

## clear the memory
rm(list = ls())
install.packages("gplots")

## Load library packages
library("plm")            # Panel data model
library("splm");          # Spatial panel data model
library("spdep");         # Define listw class, and transform it into a sparse matrix
library("lmtest")
library("gplots")
library("calibrate")      # Show labels to specific points in plots
library("xtable");        # Export results to Latex
library("texreg");        # Export results to Latex
library("orcutt")

## Import data as a data.frame.
monData <- read.csv("~/phd_thesis/Rfile/monData.csv", header = T)
wmat = as.matrix(read.csv("~/HOLA/research/green-jobs/wmat/adjmat_backup.csv", header = TRUE, row.names = 1))

monData <- read.csv("D:/research/green-jobs/data/monData.csv", header = T)
wmat = as.matrix(read.csv("D:/research/green-jobs/wmat/adjmat_backup.csv", header = TRUE, row.names = 1))

monData <- read.csv("~/monData.csv", header = T)
monData <- read.csv("~/HOLA/research/green-jobs/data/monData.csv", header = T)
qtrData <- read.csv("~/HOLA/research/green-jobs/data/qtrData.csv", header = T)

## Describe the data set
names(monData);
dim(monData);
head(monData);
tail(monData);
attach(monData);
monData$lemp <- log(emp)
monData$lwage <- log(realwage)
## monData$td <- paste(month,year,sep=",")  
write.table(monData,file="~/HOLA/research/green-jobs/data/monData.csv",,sep=",",row.names=F)

## Unit root tests
y <- data.frame(split(monData$wage,monData$county))
head(y)
purtest(y,lags = "AIC",pmax = 2, exo = "trend", test = "levinlin")
purtest(y, exo = "trend", test = "hadri")
purtest(y,lags = "Hall",pmax = 2, exo = "trend", test = "ips")
?purtest

## Create panel data frame
monPanel <- pdata.frame(monData, index = "county", drop.index = TRUE, row.names = TRUE)
names(monPanel)


## Data description
pdf("~/HOLA/research/green-jobs/latex/figures/windcap.pdf", height=6, width=8)
capmat = as.matrix(monPanel$newcap)
plot(colSums(capmat), xlab = "Month t", ylab = "Wind capacity installed", col = "red")
 dev.off()

pdf("~/HOLA/research/green-jobs/latex/figures/wells.pdf", height=6, width=8)
wellmat = as.matrix(monPanel$wells) <- data.frame(split(monData$emp, monData$county))
purtest(y,pmax = 11, exo = "intercept", test = "madwu")

sum(monPanel$wells[monPanel$year == 2011])
## Basic FDL regression model
fm.emp <-emp ~ lag(wells,0:6)+ lag(spud,0:6) + newcap
fm.emp.lrp <-  emp ~ wells + (lag(wells,1)-wells) + (lag(wells,2)-wells) + (lag(wells,3)-wells) + (lag(wells,4)-wells) + (lag(wells,5)-wells) +(lag(wells,6)-wells) 
fm.diff <-diff(emp) ~ lag(diff(wells),0:6)
fm.emp1 <-emp ~ lag(wells,0:6) + lag(wells,7:12)


## Random effects model cannot be used

job.twfe <- plm(fm.emp,monPanel,model = "within",effect = "twoways")
summary(job.twfe)
job.twre <- plm(fm.emp,monPanel,model = "random",effect = "twoways")
phtest(job.twfe,job.twre)
## ## Strict exogeneity test

## fm.diff <- diff(lemp) ~ lag(diff(wells),0) + wells 
## difftest <- plm(fm.diff,monPanel,model = "within",effect = "twoways")
## summary(difftest)

## Tests for individual and time effects 
plmtest(fm.emp,monPanel,effect = "twoways", type = "bp")
plmtest(fm.emp,monPanel,effect = "individual",type = "bp")
plmtest(fm.emp,monPanel,effect = "time",type = "bp")
## General serial correlation tests
## pbgtest(job.twfe,monPanel, order = 2)   
pdwtest(fm.emp,monPanel)
pwartest(fm.emp,monPanel,method = "arellano")
pwfdtest(fm.emp,monPanel) # test sc in differenced errors
## Cross-sectional dependence(XCD) tests (local with wmat)
pcdtest(fm.emp,monData,index = "county", w = wmat, model = "within",effect = "twoways")
pcdtest(fm.emp,monData,index = "county", model = "within",effect = "twoways")
 
## First Difference
job.fd <- plm(fm.emp,monPanel,model = "fd")
sumjobfd = summary(job.fd)
sumjobfd
job.fd1 <- plm(fm.emp1,monPanel,model = "fd")
summary(job.fd1)
## test optimal number of lags


## Robust covariance matrix estimation
vcov = vcovHC(job.fd,method = "arellano")
coeftest(job.fd, vcov)
job.gls <- pggls(fm.emp,monPanel,model = "within")
sumjobgls = summary(job.gls)
sumjobgls

LRP.fd = sum(job.fd$coefficients[2:8])/12
LRP.fd
well2011 = sum(monPanel$wells[monPanel$year == 2011])
job2011 = LRP.fd*well2011*10182150/12/6
job2011

gmmwells = c(15.77152, 13.0483251767, -2.6579219236, -9.9551518106, 5.3032277835, 8.1570373662, 11.343454618)

gmmwind = c(-0.762396, 0.0955101393, 0.0897968793, -0.1218001097, -0.1336975871, 0.0774407616, -0.3415223528)
pdf("~/phd_thesis/figures/ch3/gmmwind6.pdf", height=6, width=8)
par(xaxs="i", yaxs="i",cex.lab=1.5,cex.axis = 1.5,cex.main = 1.5)

plot(c(0:6),gmmwells, type = "l",xlab = "lag", ylim = NULL, ylab = "Job created per well completed", main = "Short term impact of shale activity on employment",font.main = 2,col = "red",lwd = 2)

lines(c(0:6),gmmwind, type = "l",xlab = "lag", ylim = NULL, ylab = "Job created per MW wind capacity", main = "Short term impact of wind activity on employment",col = "red",lwd = 2)

plot(c(0:12),job.fd$coefficients[2:8], type = "l",xlab = "lag", ylim = NULL, ylab = "Job created per well completed", main = "Short term impact of shale activity on employment",font.main = 2,col = "red",lwd = 2)

plot(c(0:6),job.fd$coefficients[9:15], type = "l",xlab = "lag", ylim = NULL , ylab = "Job created per MW wind capacity", main = "Short term impact of wind activity on employment",font.main = 2,col = "red",lwd = 2)

grid()
dev.off()

texreg(list(job.fd,job.gls))
x <- data.frame(cbind(job.fd$coefficients[2:8]*100,job.gls$coefficients[2:8]*100))
x
confd = confint(job.fd)
confd
pdf("~/phd_thesis/figures/ch3/fdwage.pdf", height=6, width=8)
par(xaxs="i", yaxs="i",cex.lab=1.5,cex.axis = 1.2,cex.main = 1.5) 
plot(c (0:6),job.fd$coefficients[9:15], type = "l",xlab = "lag", ylim = c(-0.0001,0.00001) , ylab = "Job created per MW wind capacity", main = "Short term impact of wind activity on employment",font.main = 2,col = "red",lwd = 2)
plot(c(0:6),job.fd$coefficients[2:8], type = "l",xlab = "lag", ylim = NULL, ylab = "Job created per well completed", main = "Short term impact of shale activity on employment",col = "red",lwd = 2)

plot(c(0:12),pay.fd$coefficients[2:14], type = "l",xlab = "lag", ylim = NULL , ylab = "wage increased per well/MW wind capacity", main = "Short term impact of shale/wind activity on wage",font.main = 2,col = "red",lwd = 2)
lines(c(0:12),pay.fd$coefficients[15:27], type = "l",xlab = "lag", ylim = NULL , ylab = "wage increased per MW wind capacity", main = "Short term impact of wind activity on wage",font.main = 2,col = "blue",lwd = 2)
grid()
## lines(c(0:12),confd[2:14,1],type = "l")
## lines(c(0:12),confd[2:14,2],type = "l")
## lines(c(0:6),job.fd$coefficients[9:15],lty = 2,  type = "l", col = "blue",lwd = 2)
## ## ## lines(c(0:19),SRPwcap,  type = "l", col = "red",lty = 3, lwd = 2)
## abline(h = 0, lty = 1)
## legend(4,0.0008, c("FD wells","FDGLS wells"),  lwd = c(2,2), lty = c(1,2), col = c("red","blue"))
## ## abline(v = 11, col = "red")
dev.off()

## Dynamic ADL model
fm.emp2 <-  lemp ~lag(lemp,1:2) + lag(wells,0:6)
fm.emp <- lemp ~ lag(lemp,1) + lag(wells,0:6)

## Fixed Effects
job.twfe <- plm(fm.emp6,monPanel,model = "within",effect = "twoways")
jobtwfe = summary(job.twfe)
jobtwfe
vcov = vcovHC(jobtwfe,method = "white2")
jobtwfe$coefficients[1:6]

plot(c(0:6),job.twfe$coefficients[4:10], type = "l",xlab = "lag", ylim = NULL, ylab = "Job created per well", main = "lemp ~ lag(lemp,1:3)+lag(wells,0:6)",col = "red",lwd = 2)
lines(c(0:6),job.$coefficients[3:9],lty = 2,  type = "l", col = "blue",lwd = 2)
#legend(25,8e-04, c("TWFE","FD"),  lwd = c(2,2), lty = c(1,2), col = c("red","blue"))
## serial correlation
pbgtest(fm.emp,monPanel)
pdwtest(fm.emp,monPanel)
pwartest(fm.emp,monPanel,method = "arellano")
pwfdtest(fm.emp,monPanel) # test sc in differenced errors
pwfdtest(fm.emp,monPanel, h0 = "fe") #similar to pwartest

LRP.twfe = sum(job.twfe$coefficients[3])/(1-sum(job.twfe$coefficients[1:2]))
LRP.twfe

conf.twfe

SRPwell <- c(13.58, 7.408, -7.839, -10.67792671, 9.253890384, 18.39616591, 26.82733003, 8.464103586, -18.06225939, -7.537082502, 11.47567378, 1.611168154, 0.096379617, -7.642155611, -21.18690851, -26.85203097, -10.77147597, 15.48485203, 16.30259037, 0.093188678)

SRPwcap <-  c(-0.9543, 0.121355633, 0.207863492, 0.451882154, 0.325391389, 0.077387139, -0.69609554, 0.730393886, 0.252969457, 0.098070671, 0.413418914, -0.479406086, 0.476643043, 1.091903492, 1.374381361, 1.392259515, -1.103297255, 1.824461961, 1.688131377, 2.182294983)


pdf("D:/research/green-jobs/latex/figures/twfe.pdf", height=6,
width=8) par(xaxs="i", yaxs="i") plot(c(0:19), SRP, type = "l",xlab =
"lag", xlim = NULL, ylim = NULL, ylab =
"Employment per well completed", main =
"Short term well drilling impact on employment",lwd = 2, col = "blue")
## lines(c(0:19),conf.twfe[3:22,1],lty = 3,col = "red")
## lines(c(0:19),conf.twfe[3:22,2], lty = 3, col = "red")
## segments(c(0:10),conf.twfe[3:13,1],c(0:10),conf.twfe[3:13,2], col = "red", lwd = 1,lty = 4)
## legend(2,60, c("wells", "wind capacity"),  lwd = c(1,1), lty = c(1,2), col = c("blue", "red"))
abline(h = 0, lty = 2, col = "grey")
text(2,-40,"2.5% Con. Interval")
text(6,30,"97.5% Con. Interval")
dev.off()

coeftest(job.twfe, vcov=vcovHC(job.twfe, method = "arellano", type = "HC0"))

## Graphs of fixed effects
## Spacial effects:
fe.spatial = fixef(job.twfe, type = "dmean");


pdf("~/HOLA/research/green-jobs/latex/figures/spatialEf.pdf", height=6, width=8)
par(xaxs = "i",yaxs = "i")
plot(fe.spatial,main = "County-specific fixed effects", xlab = "county #", ylab = "Fixed effects, # of employment")
textxy(fe.spatial, names(fe.spatial))
dev.off()

## Time effects:
fe.time = fixef(job.twfe, effect = "time");
par(xaxs = "i",yaxs = "i")
pdf("~/HOLA/research/green-jobs/latex/figures/timeEf.pdf", height=6, width=8)
plot(fe.time,main = "Time effects", xlab = "Month", ylab = "Time effects, , # of employment", type = 'o')
dev.off()

## Graphs
plot(job$model)
plot.new()
abline(job)
####################
## Impact on wage ##
####################
## Basic regression models
fm.wage <- lag(wage,0) ~  lag(wells,0) + lag(wells,1:12)  + lag(newcap,0)+lag(newcap,1:12)

fm.wage.lrp <- lag(realwage,0) ~ lag(cumuwells,0:11)  + lag(cumucap,0:11)
fm.wage <- realwage ~ wells + newcap + as.factor(year) + as.factor(month)
pay.fd <- plm(fm.wage,monPanel,model = "fd")
phtest(pay.fd,pay.fd0)
summary(pay.fd)
vcov = vcovHC(pay.fd,method = "arellano")
coeftest(pay.fd, vcov)
pay.gls <- pggls(fm.wage,monPanel,model = "within")
pay.gls1 <- pggls(fm.wage,monPanel,model = "fd")
pay.gls.lrp <- pggls(fm.wage.lrp,monPanel,model = "within")
summary(pay.gls)
summary(pay.gls.lrp)
LRP1 = sum(pay.fd$coefficients[2:14])
LRP2 = sum(pay.fd$coefficients[15:27])
LRP1
LRP2
texreg(pay.fd)

## Pooling effects:
pay.pool <- plm(fm.wage,monPanel,model = "pooling")
## Fixed effects:
## One-way
pay.fe <- plm(fm.wage,monPanel,model = "within")
summary(pay.fe)
pay.re <- plm(fm.wage,monPanel,model = "random")
## Two-ways
pay.twfe <- plm(fm.wage,monPanel,model = "within", effect = "twoways")
sumpay = summary(pay.twfe)
sumpay
pay.lrp <- plm(fm.wage.lrp,monPanel,model = "within", effect = "twoways")
sumpaylrp = summary(pay.lrp);
sumpaylrp
pay.twre <- plm(fm.wage,monPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(pay.twre)
waldtest(job.twfe,job.twfe1)    

pdf("~/HOLA/research/green-jobs/latex/figures/wagest1.pdf", height=6, width=8)
plot(pay.gls$coefficients[1:12], type = "l", ylim = c(-0.1,0.2), xlab = "lag", ylab = "wage increased per well/MW wind capacity installed", main = "Short term wage increased due to well drilling/wind capacity",col = "blue",)
lines(pay.gls$coefficients[13:24], type = "l", lty = 2, col = "red")
legend(2,-0.05, c("Wells", "Wind capacity"),  lwd = c(1,1), lty = c(1,2), col = c("blue", "red"))
## abline(v = 11, col = "red")
dev.off()

pdf("~/HOLA/research/green-jobs/latex/figures/wagelt.pdf", height=6, width=8)
plot(pay.lrp$coefficients[2:13], type = "l", ylim = c(-0.2,0.2), xlab ="lag", ylab = "wage increased per well/MW wind capacity installed", main = "Long term wage increased due to well/wind activity ",col = "blue",)
lines(pay.lrp$coefficients[14:25], type = "l", lty = 2, col = "red")
legend(2,-0.1, c("Wells", "Wind capacity"),  lwd = c(1,1), lty = c(1,2), col = c("blue", "red"))
dev.off()

## Tests
## Tests for individual and time effects
plmtest(pay.pool,effect = "twoways", type = "bp")
plmtest(pay.pool,effect = "individual",type = "bp")
plmtest(pay.pool,effect = "time",type = "bp")
pFtest(pay.fe,pay.pool)
pFtest(pay.twfe,pay.pool)

## Hausman test
phtest(pay.fe,pay.re)
phtest(pay.twfe,pay.twre)

## Tests of serial correlation
## joint/conditional LM test, test for ramdom effects and serial correlation under normality and homoskedasticity of the idiosyncratic errors
pbsytest(fm.wage.lag, monPanel, test = "j")
pbsytest(fm.wage.lag, monPanel, test = "re")
pbsytest(fm.wage.lag, monPanel)

## General serial correlation tests
pbgtest(fm.wage,monPanel, order = 11)
pdwtest(fm.wage,monPanel)
pwartest(fm.wage,monPanel)
pwfdtest(fm.wage,monPanel) # test sc in differenced errors
pwfdtest(fm.wage,monPanel, h0 = "fe") #similar to pwartest

## GMM  estimator

emp.gmm <- pgmm(lemp~lag(lemp,2) + lag(wells,0) + newcap | lag(lemp,13), data = monPanel, effect = "twoways", model = "twostep")
sumjobgmm = summary(emp.gmm)
sumjobgmm
lrp.gmm = sumjobgmm$coefficients[2]/sumjobgmm$coefficients[1]
lrp.gmm
emp.gmm2 <- pgmm(lemp~lag(lemp,1) + wells + cumucap | lag(emp,2:11), data = monPanel, effect = "twoways", model = "onestep", transformation = "ld")
summary(emp.gmm2)
summary(emp.gmm2,robust = TRUE)

emp.gmm <- pgmm(emp~lag(emp,1:2) + lag(cumuwells,0:6) + cumucap | lag(emp,2:132), data = monPanel, effect = "twoways", model = "twosteps")
summary(emp.gmm)

emp.gmm <- pgmm(emp~lag(emp,1:2) + cumuwells + cumucap | lag(emp,2:132), data = monPanel, effect = "twoways", model = "twosteps")
summary(emp.gmm)


#########################
## Spatial panel model ##
#########################

## ML Implementation ##
fm.emp <- emp ~ completion + lag(completion,1)
fm.wage <- lwage~wells+newcap
## Import spatial weights matrix.

## spacial weights matrix
wmatl <- mat2listw(wmat, style = "W");

## Spatial lag  model, random effects
lagre <- spml(formula = fm.emp, data = monData, index = "county", listw = wmatl,  model = "random", lag = T, spatial.error = "none")
summary(lagre)

## Spatial lag model, Two-way fixed effects
lagtwfe <- spml(formula = fm.emp, data = monData, index = "county", listw = wmatl, model = "within", lag = T, spatial.error = "none", effect = "individual", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)

summary(lagtwfe)

## Spatial error model, Two-way fixed effects
errtwfe <- spml(formula = fm.emp, data = monData, index = "county", listw = wmatl, model = "within", lag = F, spatial.error = "b", effect = "twoways",method = "eigen")

summary(errtwfe)

eff <- effects(errtwfe)

## Wage

## Spatial lag model, Random effects
lagre.wage <- spml(formula = fm.wage, data = monData, index = "county", listw = wmatl,  model = "random", lag = T, spatial.error = "none")
summary(lagre.wage)

## Spatial lag model, Two-way fixed effects
lagtwfe.wage <- spml(formula = fm.wage, data = monData, index = "county", listw = wmatl, model = "within", lag = T, spatial.error = "none", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
summary(lagtwfe.wage)

## Spatial error model, Two-way fixed effects
errtwfe.wage <- spml(formula = fm.wage, data = monData, index = "county", listw = wmatl, model = "within", lag = F, spatial.error = "b", effect = "twoways",method = "eigen")
summary(errtwfe.wage)

eff <- effects(errtwfe.wage)


#######################
## GM Implementation ##
#######################

## Spatial error model, fixed effects
GM.error <- spgm(formula = fm.emp, data = monData, index = "county", listw = wmat, moments = "initial", model = "within", spatial.error = TRUE)
summary(GM.error)

## Spatial lag model, fixed effects
GM.lagerr <- spgm(formula = fm.emp, data = monData, index = "county", lag = TRUE, listw = wmat, model = "within", spatial.error = TRUE)
summary(GM.lagerr)



#Testing#

## LM test, no random effects assuming no spatial correlation
## H0: sigma^2_mu = 0|lambda = 0
test1 <- bsktest(x = fm.emp, data = monData, index = "county",listw = wmatl, test = "LM1")
test1

## LM test, no spatial correlation assuming no random effects
## H0: lambda = 0|sigma^2_mu = 0
test2 <- bsktest(x = fm.emp, data = monData, index = "county", listw = wmatl, test = "LM2")
test2

## Joint LM test
test3 <- bsktest(x = fm.emp, data = monData,index = "county", listw = wmatl, test = "LMH")
test3

## Conditional LM_lambda test#
## H0: lambda = 0| sigma^2_mu >= 0 
testc1 = bsktest(x = fm.emp, data = monData, index = "county", listw = wmatl, test = "CLMlambda");
testc1

## Conditional LM_mu test#
## H0: sigma^2_mu = 0| lambda may or may not be zero
testc2 = bsktest(x = fm.emp, data = monData, index = "county", listw = wmatl, test = "CLMmu");
testc2
## ## Tests
## ## Tests for individual and time effects
## plmtest(job.pool,effect = "twoways", type = "bp")
## plmtest(job.pool,effect = "individual",type = "bp")
## plmtest(job.pool,effect = "time",type = "bp")
## pFtest(job.fe1,job.pool)

## ## Hausman test
## phtest(job.twfe,job.twre)

## ## Tests of serial correlation
## ## pwtest null: sigma_c^2 = 0, no unobserved effects.
## pwtest(fm.emp.fe, monPanel)
## ## joint/conditional LM test, test for ramdom effects and serial correlation under normality and homoskedasticity of the idiosyncratic errors
## pbsytest(fm.emp, monPanel, test = "j")
## pbsytest(fm.emp, monPanel, test = "re")
## pbsytest(fm.emp, monPanel)
## ## If the presence of a random effect is taken for granted
## pbltest(fm.emp, monPanel, index = "county")
## ## General serial correlation tests
## pbgtest(job.twfe,monPanel, order = 2)
## pdwtest(fm.emp,monPanel)
## pwartest(fm.emp,monPanel,method = "arellano")
## pwfdtest(fm.emp,monPanel) # test sc in differenced errors
## pwfdtest(fm.emp,monPanel, h0 = "fe") #similar to pwartest

## ## Cross-sectional dependence(XCD) tests (local with wmat)
## pcdtest(fm.emp,monData,index = "county", w = wmat, model = "within",effect = "twoways")
## pcdtest(fm.emp,monData,index = "county", model = "within",effect = "twoways")

## ## Creat identifier variable windcty/shalecty for DID method
## monData$postwind <- as.numeric(monData$cumucap > 0)
## pwmat = matrix(monData$postwind, 254,132,byrow = T)
## pwmat[Rowsums(pwmat)>0] = 1
## monData$windcty = c(t(pwmat))

## monData$postshale <- as.numeric(monData$cumuwells > 0)
## psmat = matrix(monData$postshale, 254,132,byrow = T)
## monData$shalecty = c(t(psmat))

## plotmeans(emp ~ county, main="Heterogeineity across countries", data=monPanel)

## plotmeans(emp ~ td, main="Heterogeineity across countries", data=monPanel)
